Inhomogeneous elastic response of silica glass 



F. Lconfortc, 1 A. Tanguy, 1 J. P. Wittmer, 2 and J.-L. Barrat 1 

1 Universite Lyon I Laboratoire de Physique de la Matiere Condensee et des Nanostructures; CNRS, 
UMR 5586, 43 Bvd. du 11 Nov. 1918, 69622 Villeurbanne Cedex, France 
2 Institut Charles Sadron, CNRS, 6, Rue Boussingault, 67083 Strasbourg, France 

(Dated: February 3, 2008) 

Using large scale molecular dynamics simulations we investigate the properties of the non-affine 
displacement field induced by macroscopic uniaxial deformation of amorphous silica, a strong glass 
according to AngelPs classification. We demonstrate the existence of a length scale £ characterizing 
the correlations of this field (corresponding to a volume of about 1000 atoms), and compare its 
structure to the one observed in a standard fragile model glass. The "Boson-peak" anomaly of the 
density of states can be traced back in both cases to elastic inhomogeneities on wavelengths smaller 
than £ where classical continuum elasticity becomes simply unapplicable. 
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The vibrational dynamics of glasses and in particular 
the vibrational anomaly known as the "Boson Peak" , i.e. 
an excess of the low-energy density of state in glasses 
relative to the Debye model, have attracted consider- 
able attention in condensed matter physics 0- 
This anomaly is observed in Raman and Brillouin spec- 
troscopy and inelastic neutron scattering exper- 
iments in many different systems (polymer glasses , 
silica 0, metallic glasses 0) and the corresponding ex- 
citations are often associated with heat capacity or heat 
conductivity low temperature anomalies. Many interpre- 
tations of these vibrational anomalies have been put for- 
ward, and generally involve some kind of disorder gen- 
erated inhomogeneous behavior 0, whose exact nature, 
however, is the subject of a lively debate 00, 0, Ell- 
in this work, we argue that the natural origin of these 
anomalies in "fragile" as well as "strong" glasses lies 
in the inhomogeneities of the elastic response at small 
scales, which can be characterized through the correla- 
tion length £ of the inhomogeneous or "non-ajffme" part 
of the displacement field generated in response to an elas- 
tic deformation imposed at the macroscopic scale. The 
existence of such a length has been sug gested in a sc- 
ries of previous numerical studies |lCt 111! Il2| on two and 
three dimensional Lcnnard- Jones (LJ) systems, and is 
experimentally demonstrated in macroscopic amorphous 
solids (foams [l^, emulsions 0], granulars [l^, . . . ). At 
a more microscopic level, evidence has been provided re- 
cently by UV Brillouin scattering experiments on amor- 
phous silica |l6l ] . Being a natural consequence of the dis- 
order of microscopic interactions the non-affine dis- 
placement field is responsible in particular for the break- 
down of Born-Hu ang 's formulation 0] for the prediction 
of elastic moduli |l0l Il8l [l9j , and has recently been stud- 
ied theoretically by Lemaitre et al. |l9j and DiDonna et 

al. m. 

In practice, however, it appears that the only practical 
way to quantify this effect for a given material consists in 
direct molecular simulations 0, 0] . The present con- 



tribution extends, for the first time, the numerical anal- 
ysis to a realistic model of an amorphous silca melt — 
a "strong" glass according to Angell's classification [2l| . 
Our results are compared to a previously studied "frag- 
ile" reference glass formed by weakly polydisperse LJ par- 
ticles in 3D |l2| . Strong and fragile systems have very 
different molecular organisation and bonding. Although 
the intensity of vibrational anomalies is less important 
in fragile systems, it is well documented in experiments 
on polymer glasses or in simulations of Lennard- Jones 
systems |22l . The observation of common features points 
to a universal framework for the description of low fre- 
quency vibrations in glassy systems. One recent finding 
of particular interest, is the fact that, in these LJ sys- 
tems, the Boson Peak anomaly appears to be located at 
the edge of the non-affine displacement regime, its posi- 
tion given by the pulsation associated with £ 11]. This 
begs the question whether this is a generic result applying 
also to other glasses, specifically to strong glass forming 
materials such as amorphous silica, which is character- 
ized by an intricate local packing |22( — believed widely 
to be the specific origin of the vibrational anomaly |2|. 
As in our earlier contributions, we will focus on the anal- 
ysis of the non-affine displacement field obtained in the 
linear elastic strain regime and the eigenmode density of 
states for systems at zero temperature or well below the 
glass transition. 

The amorphous silica is modelled using the force field 
proposed by van Beest et al. |23j . (For details about this 
"BKS" potential, see Refs. |24j.) We performed clas- 
sical NVT Molecular Dynamics simulations of systems 
containing N = 8016, 24048 and 42000 atoms with den- 
sity p = 2.37 g/cm 3 , in fully periodic cubic cells with 
sizes L = 48.3 A, 69.6 A and 83.8A, respectively. The 
short ranged part of the BKS potential was truncated 
and shifted at a distance of 5.5 A. For the Coulomb part 
we use the Ewald method with a real-space sum trun- 
cated and shifted at 9 A [25J. To obtain the silica glass, 
we first equilibrate all systems at T — 5000 K during 
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FIG. 1: Inhomogeneous part 5u(r) of the displacement field 
w(r) for the imposed macroscopic uniaxial strain in elonga- 
tion t xx = 5.10 -3 for a silica glass containing N = 42000 
particles (L = 83.8 A). Projection of the field is done on the 
(x — z)-plane for all particles with position r close to the plane, 
the arrow length being proportional to Su(r). The field re- 
sides in the linear elastic regime, i.e. has a magnitude varying 
linearly and reversibly with the applied deformation. As vi- 
sual inspection shows, it is strongly spatially correlated and 
involves a substantial fraction of all atoms. 

0.8 ns. An ensemble of three independent configurations 
was studied for each system size |27j. Next, we perform 
a quench from T = 5000 K to T = K by decreasing 
linearly the temperature of the external heat bath with a 
quench rate of 1.8 K/ps 24\ . Finally, a Conjugate Gradi- 
ent algorithm was used to minimize the potential energy 
of the systems yielding T = K configurations with hy- 
drostatic pressure (P) w 0.4 GPa. The static properties 
were checked against published results obtained with the 
same amorphous silica model [jjj . 

We now describe briefly the protocol used in order 
to investigate the elastic behavior at zero temperature of 
the model glasses under uniaxial deformation (for more 
details, see Refs. H El). The procedure consists in 
applying a global deformation of strain He^H < 1 to 
the sample by rescaling all coordinates in an affine man- 
ner. Starting from this affinely deformed configuration, 
the system is relaxed to the nearest energy minimum, 
keeping the shape of the simulation box constant. The 
relaxation step releases about half of the elastic energy of 
the initial affine deformation and results in the displace- 
ment Su(r) of the atoms relative to the affinely deformed 
state, defining the non-affine displacement field. A typ- 
ical field for a silica glass is presented in Fig. where 
a 2D projection of 6u(r) in the plane containing the ap- 
plied deformation direction is shown. 

This procedure allows us to measure directly the elas- 
tic coefficients from Hookc's law i- e - from the 
stress differences Aer Q/ 3 = a^f — <j™J , a r ^ being the 



total stress tensor of the reference state configuration 
(quenched stresses), and a^f the one measured in the 
deformed configuration after relaxation. From the re- 
sulting values of the Lame coefficients A — Aa yy /e xx 
and [i = (Aer^ — Ao~ yy )/2e xx one obtains the associ- 
ated transverse and longitudinal sound wave velocities, 
Ct = \/(i>/p, Cl = \/(A + 2/i)/ ( o. In the case of the sil- 
ica glass (A « 34.4 GPa, /i w 37.2 GPa, C T « 3961.4 
m/s, Cl ~ 6774.5 m/s), these quantities are in good 
agreement with data from Horbach et al. [24| and Zha 
et al. 26] (for silica under a density of 2.2g/cm 3 , taking 
into account the scaling factor (2. 37/2. 2) 1 / 2 [24| inherent 
to the choice of a higher density) . 

The linearity of the strain dependence of both the 
displacement field and the stress difference Acr Qa have 
been verified explicitly following Ref. El- The elas- 
tic (reversible) character of the applied deformation is 
checked by computing the remaining residual displace- 
ment field after removing the external strain [TlT |. An 
alternative quantification of the plastic deformation is 
obtained by considering the participation ratio Pr = 
N' 1 (Ei^) 2 /£* (Su 2 ) 2 of the noise 5u(r) El- As 
long as Pr f=a 1, all atoms are involved in the non- 
affine field, while irreversible plastic rearrangements are 
marked by Pr — > 0, with only a few particles involved. A 
choice of e xx = 10 for the L,I glass with L = 56a and 
of e xx = 5.10 -3 for the silica glass were found to ensure 
reversible and linear behavior, with 20% < Pr < 30% 
and 25% < Pr < 40%, respectively [H|. 

Visual inspection of the snapshot suggests that the 
field is strongly correlated over large distances, with the 
presence of rotational structures previously observed in 
Ref. El f° r LJ systems |2^|. In order to characterize 
this kind of structure, we normalize the field by its sec- 
ond moment, i.e. 5u(r) i— > Su(r)/ (Su(r) 2 ) 1 / 2 . In this 
way, in the linear elastic regime, it becomes independent 
of the applied strain and the system size |2(| ■ 

Next, we study the Fourier power spectrum of the 
fluctuations of this normalized field. This spectrum 
can be described by two structure factors, SL(k) = (j 
J2^ = i k-Su(r_j) exp (ik ■ r_-) || 2 ) /N relative to the longitu- 
dinal and S T (k) = (|| X)JLi kASu{ Lj ) exp (ik • r } ) \\ 2 )/N 
relative to the transverse field component |llj. These 
quantities are plotted in Fig. [21 as function of the wave- 
length A = 27r/fc, where k = kk = (2n/L)(l, m, n) with k 
being the normalized wavevector. Brackets (•) denote an 
average over the degeneracy set associated with A, and 
over an ensemble of c onfig urations. As expected from 
our study of L J glasses [llj , the longitudinal power spec- 
trum of silica is always smaller than the transverse one. 
The main difference between the two materials resides 
in the hierarchical progression of the decoupling between 
transverse and longitudinal contributions at short wave- 
lengths that appears in the case of the silica glass (the 
spectra of L J systems being only weakly wavelength de- 
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FIG. 2: (Color online) Squared amplitudes of the Fourier 
transforms Si,(k) and Sr(k) for the longitudinal (bottom) 
and transverse contributions to the normalized non-affine field 
Su(r) of silica glass at T — K under macroscopic elongation, 
plotted vs. the wavelength A = 2ir/k (in A). The various 
system sizes included demonstrate a prefect data collapse. 
The transverse contribution is more important for all wave- 
lengths. The spectra become constant for large wavelengths 
with a relative amplitude of about 10. The spectra for a LJ 
glass (spheres, length given in sphere diameters) have been 
included for comparison |3(ill 



pendent). This can be traced back to the local structure 
of silica which is represented by arrows giving the posi- 
tions of the n first neighbor shells r^_^, where n £ [1,4] 
and (a,/3) € {Si,0}. Structural effects disappear at 
distances greater than 4-5 tetrahedral units SiO^, i.e. 

r (4"5) 
{a — a} 



with a G {5*1,0}, and the longitudinal contri- 
bution to the non-affine displacement field becomes then 
about 10 times smaller than the transverse one — similar 
to our finding for LJ systems [TlLl3*oj. 

We conclude that the non-affine displacement field 
is of predominantly rotational nature in both "fragile" 
and "strong" glasses, and proceed to extract a char- 
acteristic length representative of this rotational struc- 
ture. Considering the coarse-grained field U_j(b) = 
Nj" 1 Sigy ^n{Li) °f au Nj particles contained within a 
cubic volume element V} of linear size b, we compute the 
coarse-graining function B(b) = (U_j(b) 2 }j 2 . As shown 
by Fig. we find for both glasses an exponential decay, 
well fitted by the characteristic scales £ w 23er for the 
"fragile" glass, and £ « 33 A, i.e. near 23 x rrg i _ i 
for the "strong" glass. The latter length scale has also 
been indicated in Fig. [21 The exponential behavior be- 
comes more pronounced with increasing system size (not 
shown) which reduces the regime of the cut-off observed 
at large b/L ks 1 which is expected from the symmetry 
of the total non-affine field. (B(b) — ► 1 for b — ► due to 
the normalization of the field.) 

The existence of such a characteristic length scale 
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FIG. 3: Amplitude of the coarse-graining function B(b) of 
the normalized non-affine field averaged over a volume ele- 
ment of lateral size b, versus the ratio b/L, for LJ and silica 
(squares) systems under uniaxial elongation. Since the total 
displacement is zero by symmetry B(b) must necessarily van- 
ish for large b ~ L ( "sum rule" ) . For sufficiently large system 
sizes, allowing to probe a broad a/L <C b/L -C 1 region, our 
data demonstrates an exponential decay with characteristic 
length scale £ w 23a for LJ and £ w 33 A for silica glasses. 



has already been underlined in Ref. [ll| f° r the LJ sys- 
tem, and has been related to the position of the Boson 
Peak in the density of vibrational states. In order to 
test this assumption in the case of the silica glass, we 
computed the vibrational density of states (VDOS) g(v) 
using the Fourier transform of the velocity autocorrela- 
tion function j^, calculated during 1.6 ns at T = 300 
K (followed after a run of 8 ns to assure equipartition of 
the kinetic energy at this temperature). The VDOS is 
shown in the inset of the Fig. and is in good agree- 
ment with results from [24| ■ In the main part of Fig. 0] 
reduced units x = v x £/Ct are used in order to plot the 
excess of vibrational states according to Debye's contin- 
uum prediction, i.e. g{x) / goebye{x), with £ the previous 
characteristic length scales and CV the sound velocities 
for transverse waves, for LJ and silica glasses. (The De- 
bye prediction must obviously become correct for small 
cigcnfrcqucncies. To access this frequency range even 
larger simulation boxes are needed.) This plot confirms 
the fact that the Boson Peak position can be well approx- 
imated by the frequency associated with the correlation 
length £ of elastic heterogeneities in both LJ and silica 
glasses. 

In summary, we have demonstrated the existence 
of inhomogeneous and mainly rotational rearrangements 
in the elastic response to a macroscopic deformation of 
amorphous silica. Our results are similar to the ones ob- 
tained previously for LJ glasses. The characterization of 
the non-affine displacement field demonstrates the exis- 
tence of correlated displacements of about 1000 particles 
corresponding to elastic heterogeneities of characteristic 
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FIG. 4: Inset: VDOS g{v) for the silica glass at T = 300 
K. Main figure: Excess of vibrational states g(x) compared 
to Debye's continuum model gDeb ye {x), using reduced units 
x = v£/C T , with C T ~ 4.2 (in LJ units) and C T ~ 3961.4m/s 
for LJ and silica respectively and f as indicated in the figure. 
The Boson Peak position at x ~ 1 is well approximated by 
the frequency associated with the wavelength of order £. As 
expected, the peak is more pronounced for the strong glass. 



size £ of 20 interatomic distances. The estimate of the 
frequency associated with this length is in good agree- 
ment with the Boson Peak position. The existence of 
such a characteristic length in glasses should encourage 
to view the Boson Peak as a length — rather than a fre- 
quency — marking the crossover between a regime where 
vibrations in glasses with wavelengths larger than £ can 
be well described by a classical continuum theory of elas- 
ticity, and a small wavelength regime where vibrations 
are strongly affected by elastic heterogeneities. 

In a nutshell, the vibrational anomaly is therefore sim- 
ply due to physics on scales where classical continuum 
elastic theories (such as the Debye model) must necessar- 
ily break down. This leaves unanswered the important 
question what additional excitations are probed that pro- 
duce the peak but suggests a similar description for differ- 
ent glass formers. Interestingly, the existence of a length 
scale of comparable magnitude accompanying the glass 
transition of liquids as been demonstrated very recently 
. This (dynamical) length characterizes the number of 
atoms which have to move simultaneously to allow flow 
just as our (static) length £ describes the correlated parti- 
cle displacements. Since the glass structure is essentially 
frozen at the glass transition both correlations may be 
closely related, possibly such that the non-affine displace- 
ments might be shown in future work to be reminiscent 
of the dynamical correlations at the glass transition. 

Computer time was provided by IDRIS, CINES and 
FLCHP. 
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